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ABSTRACT 


In  the  past  two  decades  since  the  advent  of  Kalman's  recursive  filter, 

■> 

numerous  algorithms  for  linear  estimation  have  emerged.  Most  of  these  algo¬ 
rithms  are  recursive  and  rely  on  solving  a  Riccati  equation  or  equivalent 
recursive  equations.  It  will  be  shown  how  some  of  the  classical  problems  such 
as  Linear  Smoothing  and  Recursive  Block  Filtering  problems  can  be  solved 
exactly  by  some  new  nonrecursive  algorithms  which  are  based  on  the  Fast  Fourier 
Transform  (FFT).  Moreover,  these  algorithms  are  readily  modified  to  generate 
the  Riccati  matrix  at  specified  times,  if  this  is  desired.  These  results  are 
then  extended  to  a  block  filtering  algorithm,  where  data  is  received  and 
smoothed  recursively  block  by  block.  Real  time  batch  processing  applications 
include  image  processing  and  array  processing  of  signals. 
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1 .  INTRODUCTION 

In  the  past  two  decades  since  the  introduction  of  Kalman's  recursive  filter 
[1,2]  several  surveys  of  this  subject  [3-5]  have  exhibited  different  algorithms 
for  the  same  basic  filter.  For  example,  Ho  and  Lee  [6]  use  a  Bayesian  approach, 
Rauch  et.  al .  [7]  utilize  a  maximum  likelihood  principle,  Meditch  [8]  uses  a 
projection  theorem,  Kailath  [9]  employs  innovations,  and  Lainiotis  [10]  uses 
the  partition  theorem  to  rederive  the  Kalman  filter. 

The  conventional  Kalman  approach  [1]  involves  the  propagation  of  a  state 
estimate  and  the  error  covariance  matrix  from  stage  to  stage.  Other  approaches 
include  the  filtering  algorithm  by  Fraser  [11]  which  is  based  on  finding  the 
information  matrix.  A  family  of  square  root  algorithms  which  recursively 
compute  the  square  roots  of  the  covariance  matrix  or  of  the  information  matrix 
is  associated  with  Potter  [12],  Dyer  and  McReynolds  [13],  Schmidt  [14],  Kaminski 
and  Bryson  [15],  and  Bierman  [16].  The  Riccati  equation  plays  a  major  part 
in  all  of  these  algorithms.  Recently,  Kailath  [17]  and  Morf  et  al.  [18], 
developed  filtering  and  smoothing  algorithms  in  which  the  Riccati  equation  is 
replaced  by  the  computationally  advantageous  Chandrasekhar  equation. 

After  Kalman  proposed  the  filtering  algorithm,  the  smoothing  problem  was 
subsequently  solved  in  the  state-space  time  domain  by  Carlton  [19],  Rauch  [20], 
Bryson  and  Frazer  [21],  Rauch  et.  al .  [17],  Meditch  [8],  Mayne  [22],  Anderson 
et.  al .  [23],  and  many  others. 


The  numerous  algorithms  that  have  been  derived  are  mostly  recursive  in 
structure  and  rely  on  solving  the  Riccati  equation  or  equivalently  the  Levinson- 
Trench  normal  equations  [24]  associated  with  the  solution  of  Toeplitz  equations. 
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In  this  report  we  present  some  new  algorithms  for  fixed  interval  smoothing, 
solution  of  Riccati  equations,  and  block  filtering  problems  that  arise  in  linear 
estimation  theory  for  discrete,  time-invariant  systems.  These  algorithms  have 
several  interesting  features  which  make  them  attractive  in  the  context  of  modern 
digital  signal  processing.  These  algorithms  are  non- recursive,  fast  and  are 
based  on  the  Fast  Fourier  Transform  (FFT).  For  example,  for  ARMA  models, 

p 

typical  recursive  algorithms  require  0(n  N)  operations,  while  the  new  algorithms 
which  utilize  the  FFT  need  only  0( (log2N+n)N)  operations  to  smooth  N+l  samples 
of  an  n  states,  single  input  single  output  ARMA  system  with  N  »  n.  Moreover, 
our  algorithms  seem  to  be  less  sensitive  to  round  off  and  truncation  errors. 

Finally  since  these  new  algorithms  utilize  the  FFT  and  are  non-recursive,  they 
could  be  used  to  process  large  data  batches  efficiently  in  parallel,  and  would 
be  well  suited  for  VLSI  architectures. 

The  smoothing  algorithm  developed  here  does  not  require  solving  the  Riccati 
equation  so  that  one  is  not  confronted  with  the  associated  numerical  problems 
such  as  insuring  the  positive  semi-definiteness  of  the  error  covariance  matrix 
e.g.,  as  in  the  square  root  algorithms.  It  is  shown  that  the  optimal  smooth 
estimate  can  be  represented  as  a  sum  of  two  components.  One  component  is  the 
output  of  a  Wiener  filter  with  discrete  frequency  response.  This  Wiener  filter 
is  associated  with  a  steady  state  periodic  system  which  is  observed  for  one 
period.  The  other  component,  called  the  boundary  response,  is  determined 
completely  by  certain  initial  and  terminal  values  of  the  observations  and  the 
Wiener  filter  output. 

Even  though  the  Riccati  matrix  is  not  required  in  our  smoothing  algorithms,  it 
plays  a  fundamental  role  in  a  large  number  of  problems.  Many  times  it  is  desired 
to  find  the  steady  state  solution  of  the  Riccati  equation  associated  with  a 
linear,  time-invariant  system.  In  some  situations,  the  given  system  is  over- 


-  3  - 


sampled  with  respect  to  the  Nyquist  rate,  e.g.,  in  radar  systems,  where  the 
signal  being  observed  may  be  a  narrow  band  signal  but  the  receiver  bandwidth 
is  much  larger.  Hence  it  may  be  desirable  to  obtain  the  Riccati  matrix  at  the 
Nyquist  rate,  i.e.,  at  equal  lags  of  time.  In  section  5,  it  is  shown  that  a 
minor  modification  of  the  smoothing  algorithm  yields  the  Riccati  matrix  at 
several  instants  of  time  via  our  FFT  approach. 

With  the  advent  of  array  processors  and  array  scanners,  data  is  often 
received  in  blocks,  batches,  packets,  arrays  or  lines.  Hence  it  is  desirable 
to  consider  filters  which  operate  sequentially  on  blocks  of  data.  Such  filter 
structures  have  been  considered  in  digital  signal  processing  for  convolution  of 
a  long  sequence  of  data  with  a  finite  impulse  response  (FIR)  filter.  We 
introduce  a  new  so  called  Recursive  Block  Filter.  As  the  name  suggests,  this 
filter  smoothes  the  data  non-recursively  within  a  block  and  recursively  from 
block  to  block.  Suppose  the  measurements  are  received  in  blocks  of  N+l  samples. 
One  is  asked  to  find  the  optimal  smoothed  estimate,  Xq,  x^1 .  ...  x^  for  i  =  1.2.. 
given  the  observations,  z^  for  k  =  0,1,  ...  iN  where  i  denotes  the  iL  block 
of  the  data  shown  below. 


It  is  an  on-line  filter  in  the  sense  that  if  one  treats  a  time-block  as 

A 

a  unit  or  a  packet  of  time,  then  one  is  asked  to  find  the  filter  estimate 
where  i  denotes  the  ith  block,  or  packet.  It  is  also  observed  that  this  filter 
is  similar  to  the  fixed  lag  smoother,  except  that  the  lag  of  this  filter  is  over 
non-overlapping  samples. 

This  recursive  block  filter  smoothes  data  block  by  block  and  can  be 
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implemented  on-line,  while  preserving  the  nice  properties  of  smoothing  as  well 
as  all  the  foregoing  advantages  of  the  new  algorithms.  Also  this  recursive 
block  filter  has  the  same  computational  complexity  as  our  new  smoothing  algorithm. 
In  fact,  by  combining  the  techniques  of  the  smoothing  algorithm  and  that  of  calcu¬ 
lating  the  Riccati  matrix  at  equal  lags,  one  easily  realizes  the  recursive 
block  filter. 

Applications  of  the  recursive  block  filter  can  be  found  in  communications 
and  telemetry  where  the  fixed  lag  smoother  is  known  to  be  useful,  in  image 
processing  where  a  block  of  data  is  available  at  one  time,  and  in  digital  on¬ 
line  deconvolution  of  finite  impulse  response  systems.  Often  in  practice,  a 
discrete  time-varying  system  is  modelled  as  a  piecewise  time  invariant  system. 

The  recursive  block  filter  can  be  extended  to  such  models  easily  by  simply 
applying  our  algorithms  to  successive  time-invariant  blocks. 

Our  results  utilize  the  fact  that  the  fixed  interval  smoother,  Kalman  filter 
and  the  associated  Riccati  equation  can  all  be  imbedded  into  a  fundamental 
boundary  value  problem.  In  Section  2  we  show  the  relationship  of  the  various 
filters  to  their  parent  boundary  value  problem.  In  Section  3  we  review  the 
permuted  controllable  canonical  form  of  state  variable  models  which  we  use  to 
develop  our  algorithms  in  Sections  4  and  5.  Section  6  contains  additional 
remarks.  The  smoothing  algorithm  for  a  special  case  is  derived  in  Section  7. 

This  provides  a  good  introduction  to  the  main  ideas  of  the  general  derivation. 

This  section  also  contains  some  numerical  examples. 
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2.  THE  FIXED  INTERVAL  SMOOTHER 
Consider  the  following  discrete,  time  invariant  system: 

(2.1)  Vl=A^  +  Bek 

(2.2)  z|<  =  Cx4(  +  nk 

where  x^e  Kn,  ekeRm,  z^cR^,  and  A,  B,  C  are  constant  matrices  of  appropriate 
dimensions. 

Assume  that  (ek)  and  (nk)  are  independent,  zero  mean,  Gaussian,  white 
processes  with  covariances 


where  6k  ^  =  0  if  k  f  Si  and  6k  k  =  1 .  Also  xQ  is  a  zero  mean  Gaussian  random 
variable  of  covariance  Qq  which  is  independent  of  Uk)  and  (nk)  . 

The  fixed  interval  smoothing  problem  for  this  model  consists  of  finding 


the  best  mean  square  estimate  of  for  k  =  0,  ...  N,  given  all  observations 
zk,  k  =  0,  ...  N. 

This  problem  is  equivalent  to  maximizing  the  conditional  probability 


p(x*lz0,  •••  ZN)  f°r  all  k. 


Applying  Bayes  rule  and  noting  that  all  random  variables  are  Gaussian,  we  get 
[e.g.,  26,27]  the  problem  of  minimizing 


,11.1  l  I  lzk  "  C**l 

'  Qn  1  k=0  *  ^ 


1  N 
1  l  I 

6  k=0 


subject  to  the  state  equations (2. 1 )  and  (2.2),  where  11*11^  denotes  the  norm 
induced  by  A.  Using  Lagrange  multipliers  this  can  be  transformed  into  the 
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following  unconstrained  problem 

N  T 

e!j"x  1  *  Jo  ^  '  BEk]  ' 

The  first  order  necessary  conditions  are  obtained  by  differentiating  with  respect 

to  V  £k  and  w 

(2.3)  Xk  -  CTR-1  (zk  -  C^)  -  ATXk+1  =  0 

(2.4)  ~  B\+1  =  0 

(2-S)  Vl  =  A*k  +  Bek 

for  k  =  1 ,  ...  N. 

Substituting  (2.4)  into  (2.5)  results  in 

(2.6)  xk+]  =  +  BKBTXk+1  ,  k  =  0 . N 

and  rearranging  (2.3)  we  find 

(2.7)  Ak  =  ATAk+1  +  cV(zk  -  C^)  ,  k  =  0 . N 

where, in  the  boundary  term  for  k  =  0 

Qq1^  -  ctr_1(z0  -  C^)  -  ata1  =  0, 

we  have  defined  Aq  -  Qq^Xq.  For  k  =  N  we  find  AN+^  =  0.  Thus  equations 
(2.6),  (2.7)  and 

(2.8)  Xq  =  QqAq  ,  An+1  =  0 

define  a  two  point  boundary  value  problem  which  is  equivalent  to  the  original 
smoothing  problem. 

A  standard  approach  to  solving  the  smoothing  problem  is  by  constructing 


the  smoothed  estimate  as  a  combination  of  a  forward  and  backward  filter  [22] 


Forward  Filter: 

(2.9)  Rkt,  ■  Ar/  *  BKBT  -  G^gJ;  R„  ■  Qq 

(2.10)  sk+]  =  Ask  +  SkHkCzk  -  Csk];  sQ  =  0 
where 

Gk  =  ARkCT  and  Hk  =  [R  +  CR)<CT]'1  . 

Backward  Filter: 

(2.11)  Xk  =  [I  '  cTHkCRk]LA\+l  +  clR_1(zk  ‘  Csk)] 

XN+1  =  0  • 

Then  the  smoothed  estimate  is  given  by 

(2.12)  £k  '  Vk  *  =k  ' 

Remarks 

i)  This  method  is  commonly  referred  to  as  the  two  sweep  method. 

ii)  Here  Rk  is  the  solution  of  the  matrix  Riccati  equation  (2.9),  and  Hk  i 
the  covariance  of  the  error  (zk  -  Csk ) . 

iii)  The  sk  obtained  from  (2.10)  is  the  one  step  predictor  E[xk+1 | ,H  <  k] 
which  arises  in  Kalman  filtering. 
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3.  REVIEW  OF  THE  PERMUTED  CONTROLLABLE  CANONICAL  FORM 

The  computations  required  in  solving  boundary  value  problems  like  (2.6)  - 
(2.7)  are  simplified  significantly  when  canonical  forms  are  used.  We  will 
concentrate  on  the  use  of  the  permuted  controllable  form  [28],  but  other  cano¬ 
nical  forms  like  the  permuted  observable  form  can  also  be  used. 

In  the  following  brief  review  it  viill  be  assumed  that  the  system 

(3.1)  x^+1  =  Ax^  +  Buk 

(3.2)  zk  =  Cxk 


is  controllable,  i.e.,  the  controllability  matrix 

(B,  AB,  ...  Anl B) 

has  maximal  rank.  The  smallest  positive  integer  y  <  n  such  that  the  matrix 


=  (B,  AB,  ...  Ay_1B) 


has  maximal  rank  is  called  the  controllability  index.  Then  the  permuted  control¬ 
lable  canonical  form  can  be  expressed  as 


(3.3) 


(3.4) 


where  we  assume  rank  B  =  m  without  loss  of  generality.  Here  the  matrices 
p..  are  projection  matrices  of  order  r^xr.^,  and  r^  <  r^.  This  can  be  viewed 
as  a  vector  ARMA  model . 

We  will  now  describe  the  change  of  coordinates  from  (3.1)  -  (3.2)  to 
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(3.3)  -  (3.4)  defined  by  the  transformations 

A  =  PAP'1  ,  B  =  PBQ"1 ,  C  =  CP'1 


and 


The  exposition  follows  [28]: 


Step  1 :  Formation  of  the  Controllability  Matrix 


(3.5)  CM  =  [B,  AB,  ...  An-1B]  . 

Step  2:  Choice  of  Linearly  Independent  Vectors 


From  CM,  a  set  of  n  independent  vectors,  forming  the  matrix  U,  are 
selected  as  follows: 


(3.6) 


U  = 


b,  b0  ...  b 
1  2  m 


Ab,  Ab,,  . . .  Ab 
i  c  m 


X  L- 

where  b..  is  the  v  column  of  B. 

Step  3:  Formation  of  State  Transformation  Matrix  P 

From  U"1 ,  a  selection  of  rows,  ei  for  i  =  1,  ...  m  is  then  made,  each 
row  corresponding  to  the  last  vector  in  each  group  of  b. .  Then  the 
transformation  matrix,  P,  is  obtained  by  multiplying  the  set  of  the 
rows  being  selected  by  A,  that  is 


As  discussed  before,  A  and  C  can  be  calculated  from 


Step  5:  Formation  of  Input  Transformation  Matrix 

The  transformation  of  input  requires  finding  the  matrix 
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Q  = 

where  each  successive  row  of  Q  corresponds  to  the  non-zero  rows  of 
the  matrix  PB. 

Step  6:  Formation  of 

Having  found  the  inverse  of  Q,  B  can  be  computed  from 

B  =  P  B  Q' 1 

Step  7:  Permutation  of  States 

The  final  form  of  the  system  is  obtained  by  permutation. 
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and 


,-l 


1  0  0  0  0 
0  0  0  1  0 
0  10  0  0 
0  0  0  0  1 
0  0  10  0 


2 

Because  A  b-j  and  A  b2  are  the  fifth  and  fourth  columns  respectively 
e^  and  e2  are  chosen  to  be  the  fifth  and  fourth  rows  of  lf^  respectively 
The  transformation  matrix,  P,  is  computed  next  as: 


and 


P  = 


e,A 


el« 


e2A 


0  0  10  0 

0  1-5  0-1 

1  -5  17  -1  2 

0  0  0  0  1 

0  0  0  1  -1 


P 


-1 


8  5  14  1 
5  10  10 
1  0  0  0  0 
0  0  0  1  1 
0  0  0  1  0 


After  the  transformation,  the  system  has  the  A,  C  matrix  as: 


= 


1 

0 

-8 

0 

0 


0 

1 

-5 

0 

0 


0 

0 

-4 

0 

0 


0 

0 

-3 

1 

-1 
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while 


0  2  1 
0  0  0 


PB 


0  0 

0  0 

1  -1 
0  0 

0  1 


A  further  transformation  of  the  input  vector  will  yield  the  desired 
system  matrix  B  as: 


B  =  PBQ 


0  0 
0  0 
1  0 
0  0 
0  1 


where 


1  1 

-ll 

-1 

1 

1 

0 

1 

0 

l 

J 

l  a 

The  final  step  is  done  as  follows:  Let  the  states  of  the  transformed 
system  be  labelled  as: 

xk  =  [*1(1)  xk(2)  xl(3)  I  5*(1)  x^(2)]T 

Therefore  by  permutation  three  subgroups  will  be  formed. 

They  are: 

xk  =  fk^  :  xk^2^  *k^ )  f  *fc(2)] 

The  system  will  have  the  form  of 
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~  1 
xk+l 


*k+l 

~  3 
xk+l 


Therefore  in  this  example 


~1 

xk 

~2 

' 

0 

1 

0 

0 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

= 

0 

0 

0 

0 

1 

xk 

+ 

0 

0 

-4 

-8 

-4 

-5 

-3 

~3 

1 

0 

0 

0 

o 

0 

1 

xk 

0 

1. 

r  i 

1 

2 

0 

1  > 

. 

L-4 

-1 

0 

0 

0J 

kk 

P]  A  [1  0]  ,  P2  A 


1  O' 
0  1 


A1  k 


^1  A 


,  a2  A 


»  c2  A 


-8  -4' 
0  0 


,  a3  A 


•  ^3  = 


5  -3] 
LO  -lj 


0 

0 


For  the  case  of  a  single  input,  single  output  system,  the  transformed 
A,  B  matrices  are  more  simple,  and  will  take  the  following  forms: 


A  = 


0  1 
0 


ala2  •••  an 


1 

0 

,  B  = 

• 

” 

0 

1 

The  method  to  find  the  transformation  matrix,  P,  in  this  case  is  slightly 
different  and  is  given  as  follows: 

Let  the  characteristic  polynomial  of  the  matrix  A  be 

A(A)  =  Det(AI-A)  =  An  +  a.|An  ^  +  ...  +  an-i^  +  an  • 

Then  from  the  controllability  matrix  [B  AB  ...  An_1B]  and  the  coefficients 
of  the  characteristic  polynomial,  the  inverse  of  the  transformation  matrix,  P”^ , 
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is  found  as 

P-'  -  [q,  q2  ...  q„] 

where 

qn  =  B 

Vl  =  Aqn  +  alqn  =  AB  +  alB 

qn-2  =  Aqn-1  +  a2qn  =  a2b  +  alAB  +  a2B 

ql  =  Aq2  +  anqn  =  A0  +  a-|An  ^B  +  ...  +  an_-|B 

Notice  that  there  is  no  need  to  perform  any  permutation  on  the  states  since 
the  transformed  system  is  already  in  the  desired  form. 


4.  THE  SOLUTION  OF  A  TWO  POINT  BOUNDARY  VALUE  PROBLEM  VIA  THE  FFT 


Now  we  turn  to  the  method  of  solving  the  boundary  value  problem  associated 
with  the  fixed  interval  smoother.  The  following  two  point  boundary  value 
problem  is  considered. 


+  cVu^cj^) 


k  =  0,  ...  N 


with  boundary  conditions 


*0  =  V^O  +  M0 

and 

^N+1  =  w 


where  Qg,  Vg  and  w  are  assumed  to  be  specified.  The  problem  (2.6)  -  (2.8)  in 
controllable  canonical  form  is  a  special  case  with  the  choices  Pg  =  0  and 
X^+1  =  0.  These  general  boundary  values  are  needed  when  the  Riccati  matrix 
and  block  filter  algorithms  are  derived.  The  are  projection  matrices  and 


C  =  (Cj,  ...  C^)  will  be  partitioned  accordingly. 

For  convenience  we  define  I  to  be  the  identity  matrix  of  order  m,  Pn  -  0, 

m  u 

Py  “  Im’  xk  =  Ak  =  Ak  =  0  for  a11  k’  ^y+1  =  -Im’  Cy+1  =  °*  and  C0  =  °* 
We  will  now  discuss  how  the  problem  (4.1)  -  (4.2)  can  be  solved  in  terms 


of  the  m  dimensional  vectors  x^  rather  than  the  n  dimensional  state  vectors 
x^.  It  follows  from  (4.1)  that 

(4.3)  x^.,  «  P.x^+1  for  i  <  y  * 


Applying  (4.3)  recursively  results  in 
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(4.4) 

where 


*k  =  Pixk-y+i  i  <  Y  ,  k  ^  y-i 


~  ^  Y_1 

Pi  =  A  Pj  ‘ 

j=i  J 

For  later  convenience,  this  expression  will  be  extended  to  all  k  >  0.  This 
can  be  done  by  using  (4.4)  as  a  definition  for  x^,  when  l  <  0.  The 
result  is  an  equivalence  class  of  vectors  for  which 

x0  =  Pix-Y+i '  1  <  Y  * 

Note  that  this  is  only  a  convenient  notation,  and  does  not  mean  that  the 
original  model  is  extended  backwards  in  time. 

Y  A 

From  (4.1)  we  get  the  1’st  component  vector  xk+1  =  xk+1  as 

i 


'k+1  "  Ai*k  +  KXk+l  * 


(4.5)  xt+1  =  \ 

i=l 

Using  (4.4),  A^+1  =  -Im,  and  =  A^P^,  where  py  =  PY  =  Im»  this  becomes 


(4.6) 


=  I  Vk-Y+i 


'k+1  "  "i^k-Y+i  +  KAk+l 


or  after  solving  for  Xk+^ : 
(4.7) 


Vi  ■  •K'1  T,  0  ‘ k  1 


Y+l 


Similarly,  equation  (4.2)  can  be  expressed  componentwise  as 

(4*8)  -k  =  AiXk+l  +  Pi-l-^k+l  +  fk  »  1  =  ]*  Y 

where  =  0  by  definition  and  fk  =  (z^-Cx^)  is  partitioned  as 

-Tn-1 


f1k  =  ciR' 


(4.9) 
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or  in  terms  of  as 


(4.10) 


fi  =  clR'1<zk  -  i,  ‘aw 


where  C.  =  C.P..  Applying  (4.8)  recursively  we  obtain  (see  (A3)) 
J  J  J 

Y-l 


(4.11) 

Substitution  of  (4.7)  and  (4.10)  into  (4.11)  results  in 

(4.12) 


Y:1  y+1  0 

J.,  ji,  *  zk  • ' 5  k  s  N-y 


where 

(4.13) 
and 

(4.14) 


°1.j  "  SjK-'Sj  +  cTr-’cj  ;cy+1 


-  YV]  rT  d*1 

Zk  '  JQ  Cy-iR  zk+i  * 


=  0 


After  an  index  transformation  (see  (Al)),  it  is  possible  to  collect  terms 
involving  x^  as 

(4.15)  J  Vk+!l  =  2k  »  1  s  k  "  N"y  • 

S.=  -Y 


Here 

(4.16) 

and 


Y+1 

I  D.  • 

j=£+l  3  K,J 


for  Z  2  0 


(4.17) 


for  Z  <  0  . 
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If  we  write  out  (4.15)  for  k  =  1,  ...  N-y  we  obtain  the  system 
~  A 


(4.18) 


1  l 

o  Yr  \ 


1 

X-y+l 

r  ~o  -| 
Z1 

• 

s 

Jo 

_XN 

_ ZN~y  _ 

By  putting  the  first  and  last  y  subvectors  on  the  right  hand  side,  we  obtain 
a  perturbed  block  Toeplitz  system: 

"  ,o 


(4.19) 


where 


'1 


’I 


L  n-y  J  L  n-y  J  L 


Ex1 

0 

0 

eV 


If  we  also  define 


■>,1 

* 

xN-y+1 

,  X  = 

- 

>  X  = 

• 

_ 

_xo  j 

_XN 

“I 

r~ 

*1 

A 

X1 

*  t  A 

xN-2y+1 

.  XY  . 

and  x  - 

Xw 

n-y 

then  the  boundary  conditions  lead  to  equations  of  the  form  (see  appendix) 

(4.20)  G^x1  =  G1^  +  z1 
and 

(4.21)  T°xt  =  T1  x1"  +  zt  . 

where  z1  and  z*  depend  on  the  observations  and  are  defined  in  the  appendix 
(also  see  Section  5).  These  equations  can  be  used  to  convert  (4.19)  into  a 
system  involving  x^ ,  ...  x^  only: 


1  1 
Z1 

E(6°)"1(G1^1  +  z1) 

0 

• 

+ 

0 

_2N-y_ 

ET(T°)'1(T1xt  +  zl) 

In  order  to  solve  this  perturbed  block  Toeplitz  system  via  the  FFT,  it  will  be 
rewritten  as  the  perturbed  block  circulant  system 


(4.23) 


Hx  =  z  +  J¥ 


where 


*my 

0  1 

1 

EIG0)"^1 

-E 

A 

J  = 

l 

O  CD 

0 

*my 

j  and  V  = 

-et 

ET(T°r1T1 

; 


j 

i 


Eqn.  (4.23)  is  solved  in  two  stages.  First  the  2mY  boundary  terms 
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are  determined,  and  then  can  be  obtained  by  circulant  deconvolution  [28], 
When  (4.23)  is  multiplied  by  we  get 


(4.24) 

b  =  JTH_12  +  jVjYb  or 

(4.25) 

b  =  [l2my  '  iVbfl-W 

Finally,  x  is  obtained  from  (4.23)  and  (4.25)  as 
(4.26)  x  =  H'h  +  H_1Jfb  . 

Eqn.  (4.25)  requires  inversion  of  only  a  2nTYx2nry  matrix;  H-1z  is  a  circular 
deconvolution  that  can  be  determined  via  the  FFT.  Once  b  is  known,  then  (4.26) 
gives  x  easily.  Since  JV b  is  a  sparse  vector  containing  only  2my  non-zero 
entries,  H-1JVb  can  be  computed  either  directly  if  2my«log2Nm  or  via  another 
FFT  based  circular  deconvolution. 

Remarks 

There  is  an  alternative  to  the  approach  presented  which  transforms 
(4.15)  into  a  perturbed  block  Toeplitz  system  by  extending  the  coefficient 
matrix  to  a  square  matrix: 


The  dimensionality  of  the  system  increases,  but  if  N  is  a  highly  composite 
number  there  may  be  advantages  in  speed  when  the  FFT  is  computed. 


1 
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5.  ALGORITHMS 

Initial ization  j 

The  proposed  algorithm  to  solve  the  two  point  boundary  value  problem  and  j 

| 

other  related  problems  will  now  be  stated.  We  start  by  summarizing  the  definitions  of 
the  necessary  quantities.  j 

i 

x,  =  x/  ,  ,  P»  =  0  ,  P^  =  I  , 

k  — k  k  — k  0  Y  m 


- — 1 
X 

1 

+ 

J 

X1 

XN-y+1 

f-  ™ 

xN-2y+l 

• 

Al  /N 

~t 

• 

-t  - 

1 

.  .  o 

X 

_ 1 

,  X  = 

X 

L  YJ 

,  X  = 

_XN 

»  X  = 

_  xn-y 

Note  that  P,.  contains  only  zeros  and  ones,  and  thus  products  with  as  a 
factor  are  not  implemented  as  matrix  multiplications. 

°i,0  :  *y*3  +  4R'1cj 


r0,£  =  ri,0  =  0  ;  ri+l ,£+l  "  ri,£  +  Di+1 ,£+l  »  0  <  i,£  s  y-1 


h  :  Y-t.v  -  •  »  s  '  s  V 


for  -y  s  i  s  0 


H  =  block  circulant  matrix  of  order  N-y  whose  first  row  is 


(Aq,  A] ,  ...  Ay,  0,  ...  0,  A_y,  ...  A_.| ) 


>  •  •  • 
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S  ^  diag(p{,  ...  P^Q"1  diag(P1  ....  P  ) 


and  let 


S  =  [diag(pj ,  ...  P^Q”1] 


Ei,* 


Gi,£  ri,£  +  Si,£ 

'  -A_Y_i+A  if  i  <  H 

0  if  i  >  a 


1  ~  T 
G  =  E 


A  A  T 

Vq  =  (0  ...  K)  ,  =  [V.A  ],  where  the  V .  are  mxn  matrices. 


Wi  =  VfCTR_1  ,  1  =  0,  ...  y-1 


T°  .  = 


<  II 

•r» 

Fi,j+1  +  WiCj  for  0  *  1 

.  ^ 

+  F-  .  if  1  <  j  <  y-i+1 

.  Fi»J 

otherwise 

' 

0 

j  <  y-i+1 

j  =  Y-i+1 
j-Y+i-1 

Ai  -Y+j 


-1+j-Y-S.Cji  otherwise 


=  E(G0)'V  ,  v12  =  -E 


f21  =  -ET  ,  =  ET(TG)-1T^ 


22 


-a 
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Find  the  inverse  of  the  block  circulant  matrix  H,  and  then  perform  the  LU- 
decomposition  of 

*  ■  'Zmy  -  jV'J’ 


(i)  Algorithm  for  the  Boundary  Value  Problem  (Figure  1) 

Step  1 :  Perform  the  FIR  filter  operation  on  the  observations  z^  to  obtain  the 


sequence 


*1  ■  %  rv-jR'\*j  •  i ;  ^  "-v 


Step  2:  Extract  the  initial  and  terminal  variables  z1  and  zb: 

k-1 


2"k  ‘  J0  Ck-jR'S  +  [~S"0]k 


t  k_I 

k  =  Wk-j-lzN-j  +  VkAfl+l  ;  1  *  k  $ 


Step  3:  Compute 


-  J 


E (G0)-1 z1 

leV)-1;* 


,  z  =  zb  +  z° 


Step  4:  Perform  the  block  circular  deconvolution 


Step  5:  Find 


l  =  H_1z 


b  = 


Step  6:  Determine  the  estimate  x  as 


x  =  y  +  H'b'Fb 


by  circular  deconvolution. 
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Step  7:  Determine  the  boundary  value  estimates 

*' .  (gW;1  * 

s'  ■  (tW^  ♦  2l] 

( i i )  The  Fixed  Interval  Smoo thincj_  Algorithm 

If  we  let  Pq  =  A^+i  =  0  in  Fig.  1,  we  arrive  at  the  fixed  interval 
smoothing  algorithm.  The  observations  are  passed  through  an  FIR  filter  with 

.  rv  _ T  _  "I 

impulse  response  C^+k  =  C  kR  ♦  -yskcO.  The  resulting  signal  (z^)  plus 
{z^}  is  deconvolved  by  H  ^  which  is  a  circular  convolution  with  the  elements 
of  the  first  column  of  the  inverse  of  the  block  circulant  matrix  H  (see 
Appendix  B).  Boundary  elements  of  this  filter  output  are  extracted  by  the 
projection  J^,  multiplied  by  boundary  filter  gain  and  then  injected  into 

a  larger  vector,  which  is  again  deconvolved  by  H_1  to  obtain  the  boundary 
response  y^.  The  final  estimate  xk  is  obtained  by  summing  the  responses  y^  and  y£. 

(iii)  Algorithm  for  Solving  the  Riccati  Equation  via  the  FFT 

The  Riccati  matrix  Rk  (see  (2.9))  can  also  be  found  from  (4.1),  (4.2), 
if  it  is  only  desired  to  compute  R^+-j  for  fixed  N.  To  derive  this  method, 
suppose  zk  =  0  for  all  k.  Then  (2.9)  implies  that  sk  =  0  for  all  k.  If  we 
set  X.  =  e.,  the  ith  unit  vector,  then  (2.12)  implies  that 

XN+1  =  RN+lei  * 


i.e.,  the  ith  column  of  R^+p  and  hence  the  matrix  R^+^  can  be  obtained  as 
solution  of  the  following  matrix  version  of  (4.1),  (4.2): 


(5.1) 


X,.,  =  AX,.  +  BKB'Ak+1  ;  Xq  =  Qf^A 


-k+1 


n~k 


0^0 


A|<  =  aTAj<+i  '  cTr  1c4  *  An+i  =  Jn 


(5.2) 
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(5.3) 


'N+l 


"  -N+l 


=  AX,+ 


BKB 


Thus,  Rn+.|  is  the  response  of  a  boundary  value  system  excited  by  an  impulse 
I  at  the  terminal  boundary.  As  a  result,  the  foregoing  algorithm  for  solving 
such  boundary  value  problems  can  also  be  used  to  construct  the  Riccati  matrix 
R^  columnwise. 

Conceptually,  it  appears  that  n  consecutive  vector  boundary  value  problems 
have  to  be  solved.  However,  there  are  some  significant  simplifications  which 
reduce  the  problem  to  solving  one  circular  deconvolution  problem  of  size  0(N) 
and  one  matrix  inversion  of  size  2my.  The  steps  1  through  7  given  above  are 
replaced  by  the  following  steps  obtained  by  setting  z^  =  0  in  the  main  algorithm. 
In  the  following  we  will  use  uppercase  variables  to  signify  the  fact  that  X,  A, 
etc.  are  now  matrices. 

Step  1 :  Extract  the  terminal  variables 


Step  2: 


Step  3: 


Determine  the  terminal  y  submatrices  Z*  of  the  block  vector  Z 

2l  =  ET(T°)'^Zt  . 


Compute  the  block  vector  J  Y  as 


JTY  *  jVb 


Zl 


Note  that  only  the  boundary  values  are  needed  in  the  sequel,  hence  direct  con¬ 
volution  may  be  more  efficient  than  circulant  deconvolution  via  the  FFT.  The 
order  of  J^fT^J  is  (2nry)  x  (2rrry). 

Step  4:  Find  the  boundary  values 

B  =  O'Vy 

where  <J>  =  C 1 2rrry  '  as  before. 


-  28 


,  so  the  last  components  are  obtained  from 

-t  =  (T0)-lT1xt  +  (tW  (see  (4.21)) 

Step  6:  Using  (4.21)  and  the  elements  of  we  can  obtain  X^»  the  state 

estimate  at  time  k  =  N.  This  gives  R^  via  (5.3).  This  algorithm  is  illustrated 

in  Figure  2.  Steps  1,  2  and  3  above  are  equivalent  to  circular  convolution 

of  the  matrix  sequence  [0,  0,  ...,  0,  Z^,  Z^,  ....  Z^]  with  [//  ^]k  and  extract’  :i 

only  the  boundary  values  of  the  output  Y^.  The  FFT  algorithm  is  needed  only 

once  to  compute  [H*^  1^,  which  are  the  matrix  block  elements  of  the  first  column 

of  the  block  circulant  matrix  H  \  Assuming  that  the  systems  in  steps  4  and  5 

are  solved  by  LU  decomposition,  and  the  LU  decompositions  of  f  and  T°  as  well 
1  T  -1 

as  the  matrices  T  and  J  H  J  are  precomputed,  the  computational  effort  for 

3 

obtaining  will  be  approximately  6n  ,  and  does  not  depend  on  the  number 

N  which  compares  favorably  with  the  direct  recursive  method  (2.9)  -  (2.10) 

o 

which  has  a  complexity  of  0(Nn  )  when  A  is  in  the  given  canonical  form.  The 
savings  are  even  more  dramatic,  when  the  algorithm  is  used  to  compute 
Rp(N+l),  i  =  1,  ...  L,  since  only  the  matrices  V  and  4>  will  have  to  be  updated 
at  each  step.  We  will  elaborate  on  this  point  when  we  discuss  a  recursive  block 
filter  in  the  next  paragraph. 

(i v)  Recursive  Block  Filtering 

As  an  application  of  the  ideas  presented  so  far,  we  turn  to  a  recursive 
block  filter,  i.e.,  a  filter  which  smoothes  one  block  of  data  at  a  time. 

This  will  be  an  approximation  to  the  global  smoother,  and  can  be  used 
more  readily,  when  data  is  only  available  in  blocks,  or  when  a  time  variant 
system  is  modeled  by  piecewise  time  invariant  systems. 

The  block  smoother  operates  just  like  the  global  smoother,  except  for  the 
initial  conditions.  If  the  blocks  are  indexed  by  r,  then  the  initial  covariance 
for  the  block  r  is  Rn+i  r_1 ,  assuming  the  block  size  is  N,  where  R^+1  r_^  denotes 


Step  5:  Recall  that  B  = 
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the  Riccati  matrix.  The  initial  value  Xg  r  also  has  to  include  the  one  step 
predictor  s^  r_-|  from  the  block  (r-1).  Thus,  the  block  filter  is  of  the  form 


(5.1) 

(5.2) 

(5.3) 


4+1  ,r  =  A4,r  +  BKB  4+1  ,r 

W  ‘  AVl,r  +  ‘  C*k,r’ 


4,r  =  Q0,r4,r  +  ur  ;  4+1, r  =  0  • 


Q0,r  =  RN+1 ,r-l ‘  yr  =  sN+l,r-l 


This  algorithm  uses  the  two  algorithms  previously  described  recursively. 

Step  1 :  Set  r  =  1,  initialize  the  parameters  by  setting  Qg  ^  =  Qg,  s^+^  g  =  0. 
For  block  r,  process  the  data  using  the  boundary  value  algorithm  for  (5.1)-(S.3). 
Step  2:  Find  the  solution  RN+-|  r  of  the  Riccati  equation  at  the  end  of  the 
data  block  using  algorithm  (iii). 

Step  3:  Update  the  following  matrices: 

Q0,r+1  =  RN+l,r 

ur+l  =  SN+1  ,r  =  A4,r 
then  use  Qg  r+^  to  update  GB,  0,  and  4'. 

Step  4:  If  all  blocks  are  processed,  stop,  otherwise  set  r  =  r+1  and  go  to 
Step  2. 

This  algorithm  is  illustrated  in  Figure  3. 
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REMARKS  AND  COMPUTATIONAL  COMPLEXITY 


a)  We  note  that  the  above  smoothing  algorithm  is  nonrecursive  and  does  not 
require  the  solution  of  a  Riccati  equation  or  equivalent  recursive  computa¬ 
tions  such  as  the  solution  of  Toeplitz  systems  via  Levinson-Trench  type 
algorithms  [24,25]. 

b)  From  Fig.  1  the  smoothing  filter  output  can  be  written  as 


(6.1) 

(6.2) 


X  =  y  +  vb  =  X°  +  Xb 

xk  yk  yk  xk  xk 


^  LH_1zb3k  +  y^ 


The  solution  component  { x°)  could  be  considered  as  the  FFT  Wiener  smoothing 
filter  output,  which  is  obtained  by  first  sampling  the  Wiener  smoothing 
filter  output  in  the  frequency  domain  followed  by  its  inverse  Fourier 
transform.  The  frequency  domain  Wiener  filter  equation  is  obtained  by 
considering  an  infinite  duration  filter  in  the  steady  state.  Specifically, 
the  steady  state,  dynamic  system  equation  can  be  written  as 


(6.3) 


Jtxk-Y+«,-l  ”  ek-l 


(M)  2k  ■  j,  Wy*  *  "k 

where  kc(-00,00).  and 

od 

(6.5)  X(u>)  =  l  xkexp(-jku)) 

k=-°°  K 

denotes  the  Fourier  transform  of  (xk).  Then  the  frequency  domain  Wiener 
filter  estimate  X°(u>)  for  the  smooth  estimate  of  is 

(6.6)  X°(u>)  =  CC*(co)R_1  C(co)  +  S^1  (a))]'1C*(io)R'1Z((o) 
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where  Sx(u>)  and  R  denote  the  power  spectrum  densities  of  xk  and  nk 
respectively,  *  denotes  the  conjugate  transpose,  and 

(6.7)  C(w)  =  \  C{exp[-j(Y-£)to] 

£=1  * 

The  quantity  S~^ (w)  is  obtained  from  the  state  equations  (6.3),  (6.4)  as 

Y+l 

(6.8)  S"1^)  =  l  l  AlK_1A.exp[j(k^)w] 

x  k=l  «.=  1  *  K 

Now  defining 

(6.9)  H(w)  =  C*(o))R'1C(o))  +  S"1  (a») ,  Z°(co)  =  C*(w) R' 1  Z(o>) 


we  can  write  (6.6)  as 

(6.10)  H(<d)X°(o>)  =  Z°(a>) 

Ort 

which  gives  the  discrete  Fourier  transform  of  (4.12)  at  u  *  , 

n  =  0,  . . . ,  N-y-1 . 

2  nn 

Thus  if  we  sample  (6.6)  at  w  =  -jp"  ,  n  =  0,  ...  M-l ,  i.e.,  let 

M-l  „ 

ZM  .  I  i  S(»  -  ?E) 
n=0 

*V)  =  Y  x°«(w  -  ^f) 

n=0  "  •' 

in  (6.6),  and  take  the  inverse  Fourier  transform,  we  obtain  the  solution 
of  the  circular  deconvolution  problem  defined  by 

x°  =  H  ^2  for  M  =  N-y 


34  - 


Moreover,  it  is  easy  to  verify  that  {x0}  and  {z  }  are  the  M  step  inverse  discrete 

n  n 

Fourier  transforms  of  { x1 2^ >  and  {z^}  respectively,  and  could  be  implemented 
via  FFT. 

c)  If  the  system  is  not  controllable,  but  remains  observable,  the  above  algorithm 
can  be  modified  by  working  with  a  canonical  observable  model.  Details  are 
left  to  the  reader. 

d)  The  algorithm  given  above  is  based  on  an  extension  of  a  method  of  inversion 
of  banded  Toeplitz  matrices  to  banded  block-Toeplitz-matrices  [29 J. 

e)  In  the  special  case  of  single  input  systems,  the  dimensionality  is  signifi¬ 
cantly  reduced.  Moreover,  the  block  circulant  matrix  H  becomes  a  symmetric 
circulant  matrix,  and  thus  also  the  Fourier  transform  of  the  first  row 

of  H  will  be  real,  which  reduces  the  number  of  operations. 


Computational  Complexity 

We  discuss  separately  the  initial  computational  effort,  which  can  be 
computed  off-line.  Secondly,  we  address  the  complexity  of  the  data  processing. 
Assume  the  input  dimension  is  m,  the  state  dimension  is  n,  the  output  dimension 
is  p. 

Initial  effort:  (let  M  =  N-y) 


(1 )  Find  K"1  and  R~^  in 

(2)  Find  K_1A.,  R_1C.  in 

J  J 

(3)  Ajlf^,  cJr_1C. 

(4)  y  matrices  IT  are  needed 

(5)  F  requires 

(6)  Set  up  f  in 


m^  +  p^  op. 

2  2 

n(m  +  p  )  ops. 

(p3 4 5 6  *  ops. 

2yir  ops. 

2-ym3  ops. 

fllSY)3  *  ♦  V3  -  J1 
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(7) 

(8) 


Find  the  LU  decomposition  of  [I  - 


( 2rrry )  2my 
3  ‘  3 


Find  H' 


2m2MlogM  +  Mm3 


(9)  _ 8y3m3  operations _ 

2  3 

Total :  M(2m2logM  +  m3)  +  m3(l  +  |^y3  +  ^  ^-)  +  m2n  +  np2  +  P  ^  -  my 


2 

Nnp  +  y  mp  ops. 
o 

4(my)  ops. 

? 

2mMlogM  +  Mm  ops. 

2  •  (2my)2 

2mMlogM  +  Mn2  ops. 

2 

4(my)  ops. 

Nnp  +  4mMlogM+2Mm2  +  16m2y2 


These  operations  are  for  the  controllable  canonical  form.  Additional  operations 
are  required  for  transformation  of  the  system  parameters  and  state  estimates 
if  the  given  system  is  not  in  this  form. 
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7.  EXAMPLES 

Example  1:  To  illustrate  the  general  ideas  of  the  derivation  of  the  smoothing 
algorithm  via  FFT  we  consider  the  special  case  of  (4.1)  -  (4.2)  where  m  =  n, 

Y  =  1  and  the  covariances  K,  Qq,  and  R  are  identity  matrices. 

(7.D  xk+1  =  Axk  +  Xk+1 

(7.2)  Ak  =  ATAk+]  +  CT(zk-Cxk),  k  =  0,  ...  N 
with  boundary  conditions 

x0  =  X0  and  XN+1  =  0  ’ 

When  (7.1)  is  solved  for  Xk+^  and  the  expression  is  substituted  into  (7.2) 
we  obtain 

(7.3)  xk  -  Axk_-|  -  A  C x i  —Ax k )  +  C  (zk~Cxk),  k  -  1,  ...  N-l  . 

We  can  collect  terms  to  get 

(7.4)  -Ax|<_1  +  (I+ATA+CTC)xk  -  ATxk+]  =  CTzk  . 


If  we  define 

A_]  =  -A  ,  AQ  -  (I+ATA+CTC)  ,  =  -AT 


then  we  can  write  the  resulting  system  (7.4)  for  k  =  1,  ...  N-l  as 


(7.5) 


By  putting  the  boundary  terms  involving  Xg  and  x^  on  the  right  hand  side  we 
obtain  the  perturbed  block  Toeplitz  system 


*'*"'#**&  4*+ jflr  V 
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(7.6) 


Now  x0  and  can  be  eliminated  by  applying  the  boundary  conditions. 

The  initial  condition  Xg  =  Ag  combined  with  (7.2)  leads  to 

(7.7)  Xg  =  ATA^  +  CT(Zg-CXg)  . 

(7.1)  can  be  used  to  eliminate  A-| . 

(7.8)  Xg  =  AT(x^-Axq)  +  CT(Zg-Cx0). 

Collecting  terms  we  obtain  the  initial  relationship 

(7.9)  (I+ATA+CTC)xq  =  ATx1  +  CTz0 

which  is  clearly  uniquely  solvable  for  xQ  in  terms  of  x1  and  zQ. 

The  terminal  condition  A^  =  0  combined  with  (7.2)  leads  to 

(7.10)  \N  =  CT(zn-Cxn)  . 

(7.1)  can  be  used  to  eliminate  A^: 

(7.11)  x^  -  Ax^j_-j  -  C  (z^-Cx^)  . 

Collecting  terms  we  obtain 

(7.12)  (I+CTC)xn  =  AxN-1  +  CTzn 

which  is  again  uniquely  solvable  for  x^  in  terms  of  x^_^  and  z^ 

Thus  (7.9)  and  (7.12)  can  be  used  to  eliminate  Xg  and  x^  in  (7.6):  let 
G0  =  I+ATA+CTC  and  TQ  ^  I+CTC,  then 
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(7.13) 


A-k 

A.  0 

S> 

r«,  ] 

» 

4 

f 

O  1—  %  • 

1 _ 

+ 

-A_16‘1(ATx1+CTz0) 

0 

• 

0 

0 

-1  0_ 

_XN-1  _ 

~o 

uZN-l_ 

_-Vo  <Axn-i+c  zh]  _ 

Finally,  this  problem  can  be  rewritten  as  a  perturbed  block  circulant  system 


(7.14) 


=  z  +  Jy 


,  where 


LXN-lJ 

L 

XN-1  _ 

X 

o 

/ 

><£ 

0 

A.r 

A  i 

0 

„  0 
_A1 

\\ 

0 

0 

\ 

-1 

si1 

V 

,  1  = 

2°  + 

* 

0 

_-¥o'c\  _ 

I 

0 

A,  _ 

m 

0 

0 

,  T  = 

-A,To’A_ 

0 

I 

?  <T" 

_ 1 

— 

letting  b  = 


we  can  first  solve  for  the  boundary  term  b  as 


(7.15)  b  =  jV1!  +  jVXb  . 

This  is  a  system  of  order  2m.  Finally  the  components  x^  ...  xN_2  can  be 
obtained  via  FFT  as 


=  H_1  (z+J'Fb)  . 
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Example  2  (A  detailed  solution  of  Riccati  Equation) 


Suppose  the  Riccati  matrix  Rg  has  to  be  evaluated  for  the  following  systems 
parameters: 


_  0 

1 

0  " 

‘  1  0" 

A  = 

-1 

-2.1 

-1.2 

-2.2 

-1.3 

-2.3 

K  = 

_0  1  _ 

C  = 


[' 


1  0  0 
0  1  0 
0  0  1 


i.e.,  this  is  a  system  with  two  inputs  and  one  output.  Note  y  =  2  here. 
We  obtain 


A 


0 


13.69  6.62 
6.62  8.98 


7.02  9.33 
1.3  2.30 


"1  2.1" 
0  0 


7.41 

0 

5.82 

7.13“ 

0 

1 

0 

0 

5.82 

0 

13.69 

6.62 

7.13 

0 

6.62 

8.98 

-1 

-2.1 

0  0“ 

0 

0 

0  0 

-7.02 

-9.33 

-1  -2.1 

-1.3 

-2.3 

0  0 

1.2 

1.3  1 

°  1 

3.2 

2.3  0 

2 

-  .2 

0  0 

-1.2 

-2.3 

2  0 

-2.3 
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"0  0-1  0 
1  0  0  -2.1  0 
-1  0  -1.2  -1.3 

_ -2. 1  0  -3.2  -2.3 


The  first  block  row  of  the  block  circulant  matrix  is  obtained  as 

w-l=[~.594  -.181  -.124  -.479  .059  .081  .059  -.077  -.124  -.OOf 

1»*  [-.181  .656  -.003  .134  -.077  .045  .081  .045  -.479  .134. 


PHI 


PSI 


.983 

-.008 

-.004 

-.009 

.071 

0.000 

-.202 

-.031 

.024 

1.029 

.001 

.002 

-.119 

0.000 

-.752 

-.051 

.050 

.041 

1.007 

.014 

-.000 

0.000 

.152 

.022 

-.008 

-.009 

.001 

1 .001 

.001 

0.000 

.142 

.089 

.256 

.188 

.024 

.051 

1.070 

0.000 

.105 

.091 

-.003 

.058 

-.002 

-.004 

.020 

1.000 

.058 

.035 

-.714 

-.897 

-.114 

-.239 

-.160 

0.000 

.811 

-.199 

.143 

.219 

.064 

.134 

-.069 

0.000 

-.093 

.912 

4.462 

5.639 

.704 

1.478 

1.000 

0.000 

7. 

.020 

1. 

300 

5.639 

7.232 

.875 

1.837 

2.100 

0.000 

9. 

,330 

2. 

300 

.704 

.875 

.116 

.243 

0.000 

0.000 

1. 

,000 

0. 

000 

1.478 

1.837 

.243 

.511 

0.000 

0.000 

2. 

,100 

0. 

000 

1.000 

2.100 

0.000 

0. 000 

1.057 

0.000 

2. 

.766 

1. 

217 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

0.000 

7.020 

9.330 

1.000 

2.100 

2.766 

0.000 

10.471 

3.153 

1.300 

2.300 

0.000 

0.000 

1.217 

0.000 

3.153 

1.429 
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"o  i  o  j  ro  o 

0  1  1  ,  B  =  1  0 

12  -5  12  0  1 


C  =  [l  0  1  ]  , 

K  =  f  1  0  1  ,  R  =  .16,  Qn  =  0 

LO  .49  J  °  |_0 


.05  0  0 

0  .01  0 
0  0  .09 


This  is  an  unstable  dual  input  single  output  system, 
a)  When  the  FFT  method  is  applied  with  N  =  7  we  obtain 


2.3035  1.5797  -10.807 

1.5797  2.7446  -  6.2513 

-10.807  -6.2513  70.4289 


R16*  R8’  i-e-’  this  is  approximately  the  steady  state  solution, 
b)  The  two  sweep  method  gives 


1.988  1.307  -9.250 
1.307  2.248  -4.667 
-9.250  -4.667  67.520 


2.166 

1.494 

-10.156 

1.492 

2.677 

-  5.825 

-10.175 

-5.821 

67.380 

2.018 

1.531  - 

9.788  ' 

6.157 

3.538  - 

24.617 

41.746 

3.024  - 

139.93 

contains  elements  of  order  10°,  which  shows  the  divergence  of  this 
method  due  to  instability. 


Figure  4:  Riccati  Equation  Solution  via  the  FFT  Method, 
This  example  shows  the  numerical  stability  of 
the  FFT  algorithm. 
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Figure  4  shows  the  evolution  of  the  element  (1,1).  The  equally  spaced 
dashes  show  the  sampled  values  using  the  FFT  approach,  and  the  other  line 
shows  the  unstable  behavior  of  the  recursive  method.  Note  that  the  line  had 
to  be  clipped. 

Example  4:  (Smoothing  and  Recursive  Block  Filtering) 

The  following  two  examples  illustrate  solutions  to  the  linear  smoothing  and 
recursive  block  filter  problems.  The  performance  of  the  filter  is  compared  for 
different  block  sizes.  The  following  system  parameters  were  chosen: 


Example  4.1 :  m=l,n=2,p=2 


(see  figures  5-8). 


Example  4.2:  m=l,n=2,p=2 


f  0  1] 

P  = 

M  •51  o  = 

.36  O' 

1 

ro 

[.5  lj»  Q0 

0  .49 

0 

1 


(see  figures  9-12). 

Note  that  the  variances  of  the  noises  are  considerably  higher  in  Example  4.2. 

For  both  examples  the  results  are  presented  in  the  same  format.  The  solid 
curves  display  the  exact  values  of  the  states  xk  =  x(<(2),  and  the  dotted  curves 
refer  to  the  estimates.  Thirty  data  points  are  shown  in  each  case.  The 
following  table  contains  the  pertinent  information.  MSE  refers  to  the  mean 
square  error  between  the  exact  and  estimated  values.  Figures  7,  8  and  10,  11 
are  recursive  block  filters  with  different  block  sizes.  Notice  that  the  block 
filter  estimates  are  quite  close  to  the  larger  interval  smoothing  filter 
estimates. 


Example  4.1 


Figure 

number  of 
observations  used 

number  of 
blocks 

block  size 

MSE 

5 

51 

1 

51 

.4497 

6 

30 

1 

30 

.7448 

_ 

30 

2 

15 

1.60 

8 

30 

3 

10 

1.099 

Figure 

number  of 
observations  used 

number  of 
bl ocks 

block  size 

MSE 

9 

51 

1 

51 

3.711 

10 

30 

1 

30 

3.811 

11 

30 

2 

15 

4.872 

12 

30 

3 

10 

4.115 

Example  5:  (Sampling  of  a  Riccati  Matrix) 

This  example  shows  the  evolution  of  the  Riccati  matrix  for  the  following 
systems  parameters:  m=l,n=2,  p=l 


A  = 


(-.9?  -1 .94) ’  ^0  =  [l  2).  C-  0.1).  B-  (°  ).  R-  .2.  K-  .1. 


Figures  13-15  display  the  components  of  Rk  using  the  two  sweep  method  as  well 
as  the  sampled  values  using  the  FFT  method.  Note  that  the  FFT  method  gives  the 
same  values  as  the  recursive  method  at  the  sampled  instants. 
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Conclusions 

Fast  nonrecursive  algorithms  for  solving  the  discrete  time  fixed  interval 
smoothing  problem  have  been  presented  along  with  extensions  to  determine 
samples  of  the  associated  Riccati  matrices  at  fixed  times.  Both  algorithms 
have  proven  to  be  useful  in  implementing  a  block  recursive  filter  to  process 
data  in  batches. 

The  advantage  of  these  algorithms  does  not  merely  lie  in  a  reduction  of 
the  computational  complexity  compared  to  more  conventional  approaches.  Their 
nonrecursive  structure  allows  parallel  architectures  and  the  use  of  FFT 
firmware  which  results  in  faster  execution  times,  and  possibly  a  reduction 
in  the  accumulated  roundoff  error. 


Figure  6:  Fixed  Interval  Smoothing,  block  size  30,  Example  4.1. 


Actual  States 


Figure  8:  Recursive  Block  Filter,  block  size  10,  Example  4.1 


Figure  10:  Fixed  Interval  Smoothing,  block  size  30,  Example  4.2 


igure  12:  Recursive 


FFT  Method 


Figure  13:  Riccati  Equation  Solutions 


Figure  14:  Riccati  Equation  Solutions 


FFT  Method 


Figure  15:  Riccati  Equation  Solution 


APPENDIX  A 


A.l  A  Technical  Lemma: 

Lemma  (Change  of  Indices  of  Summation  over  a  Rectangular  Grid): 


.  J.«ij 
'l  1=11 


i2+jrl  £-i- 

I  l 


W1 


1  L  a£,_i  i  +  1  L  a£_i  i  +  L  1  £ 


provided  >  ( ’ 2-i  1 ) *  w^ere  ^  ~  i+J. 

Proof  For  simplicity  let  i_  =  i - i -j  and  j_  =  j-j.. 

Then  by  assumption  i2  >  Let  £  =  i_+jk 

(See  Figure  16). 

Then 

for  £  <  i_2  :  0  <  j_  <  Z 

for  i^  <  £  <  j^  :  £-l2  <  i  <  l 

for  £  i  i2  :  i-lg  <  i  <  i2  Figure  16 

describes  the  range  for  j_.  Note  that  in  the  special  case  i^  =  i2  the  middle 
formula  is  redundant.  We  now  have  to  find  the  corresponding  limits  in  terms 
of  i , j  and  £  =  i+j.  Since  £  =  £+i-]+j-|»  or  £  =  £-i-|-j^  these  conditions  become 

For  *  i2-i'1  :  0  s  j-j1  s  £-ij-jj 

For  12-t,  <  s  j2-J,  :  t-vV'W 4  4  ‘"VA 

For  2  j2-j,  :  t-1,-jr02-t,)  s  J-J,  5  J2-j, 

This  can  be  simplified  to: 

For  £  <  j  1  +i 2  :  j-j  <  j  <  £-1^ 

For  i2+j1  s  £  <  j 2+i  1  :  £-i2  s  j  <  £-i-j 

For  £  >  j2+i1  :  £-i2  <  j  <  j2 

Again,  if  i 2 - i ^  =  J2~Ji »  the  middle  formula  is  redundant.  This  proves  the 
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result  and  also  its  following  corrollary: 
Corrollary :  if  l^-i-j  =  j 2-^ 1  ’  then 

j2  i2  i]+j2  £-i1 


(A2) 


l  l  a 

=J1  i=i 


I  .  I.  a£_j  j 


i2+j2 

I 

it=  i  i  +  j  2+1 


j2 


l 

j=t-i 


aP--j ,  j 

2 


A. 2  Initial 
Lenina : 

(A3) 


Conditions : 

4  =  AlXk+l 
4  =  AiXk+l 


+  f, 


i-I 

+  f{+  1 
k  j-1 


i-1 
n  i 
£=i-j 


V-jVk+i +  4;j) 


Proof:  The  case  i=l  follows  from  (4.2).  Then  use  induction.  For  i=2  we 
obtain  from  (4.2): 


4  =  A2Xk+l  +  44+1  +  fk 
=  A2Ak+l  +  PTATXk+2  +  44+1  +  fk  » 

which  is  the  claim  for  i=2.  Now  assume  the  result  is  correct  for  i.  Then 
from  (4.2)  we  find 


,i+l  -  nl*\ i  .  flT  ,  ,  -ri+1 

4  Hi4+1  Mi+lAk+l  Tk 


By  the  induction  hypothesis  we  can  replace  .  It  suffices  to  show  that 


rV 

Ki4+1 


i 

I 

j  =  l 


n 


£=i- j+1 


|(Ai-j+lXk+j+l  +  fk+j+1) 


Now 
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i-1 


Pi^k+1  =  Pi^Aixk+2  +  fk+l  +  X 

J  ■ 


i  “1  y 

n  pJ 

2.=  i-j 


jAj+k+2  *  fk»j*l'] 


T  T  T  i  ^ "1 

=  PiAlx,^  +  pjfi.,  +  l 


i  i  k+2  i  k+1 


j=1 


1  y 

n  p] 

£.=  i-j 


(Ai-jXk+j+2  +  fk+j-*X 


Pi (AIXk+2  +  fk+l }  +  -l2 


n  p 

£=i- j+1 


,T  x  +  fi-j+l 
M-j+lAj+k+l  Tk+j 


i 

i 

j=l 


n  P, 


[n=i+l-j 
which  had  to  be  shown. 


aT  ,  + 

Ai+1 -j  j+k+1  Tk+j 


i-1 


i  i  » 

With  the  convention  n  P.  -  I,  we  can  write  (A3)  for  k=0  as 

2=i 


i-1 


Aj=  l 

^  j=0 


y 

n  pi 

£=  i  -  j 


*T-jVi *  fP  • 


Now  we  can  use  (4.7)  and  (4.10)  to  obtain 


4 =  % 


i-1  1 

n  p' 

2=  i  -  j  * 


{AI-j('K  1  J1Xxj-Y+k)  +  Cl-jR  1(zj  '  J/kVj-Y+k)}- 


Using  (4.13)  and  (4.14)  this  gives 


n  pTxI  =  l  l  -  0.  .  kx.  .  +  l  Cj.1R“'zj. 

£=i  ^  j=0  k=l  J,Kn  j=0  1  J  J 


i-1  y+1 


i-1 


t-T  b-1. 


In  order  to  collect  terms  in  we  have  to  change  indices.  Use  (A1 )  and 
the  convention 


H=Y+1 


a,.  =  0 
2.1 


to  arrive  at  (j^  =  1,  j2  =  Y+1 »  iy  =  °»  ^  =  i-l) 


-  A4  - 


iipT^  •  i,  j,  •  * 


Y  A 

+  l  l  -  D 

£=i+l  k=£+l-i 


i-wk.n-y 


(A4) 


Y+i  y+1 

+  I  I  ^i-£+k  kx£-> 

£=y+1  k=£+l-i  1  *  K,K  *  ^ 


On  the  other  hand,  from  the  initial  condition  we  know  that 


Xq  “  Qo*o  +  % 


or  using  (4.4) 


^0  "  V 


Plx-Y+1 


PY*0 


*  Qo\ 


yQ  =  0  for  the  smoothing  problem.) 

S  =  diag(p|...P^)QQ1diag(P1...PY), 

S  =  [diag(p|. . .P^)Qg^ ] 

These  definitions  allow  us  to  write  out  equations  for  x_y+^...Xq,  which  in 
turn  yield  Xq  uniquely.  Then  we  can  combine  this  with  (A4)  to  obtain  a  system 
of  the  form 


(Note: 

Let 

where 


(A5) 


g0;1 


-  g1;1 


+  z 


where 


-  A6  - 


A. 3  Terminal  Condition 


Before  writing  out  the  terminal  equations,  some  preliminaries  are  needed. 
Lemma: 


(A10) 


A*  - 


-i =  +  (AT)i+Vi 


Proof:  The  case  i=0  follows  from  (4.2)  when  k=N.  Suppose  the  result  holds 
for  i.  Then  (4.2)  and  the  induction  hypothesis  result  in 


^N-(i+l)  A  ^N-i  +  fN-i-l 


'  +  <aT>1+%W  +  fN-i-l 

vJ  ^ 

■  2(aT>h,'\-  *  (aT>i+v, 


Let 


(All) 


A  = 


r°<l\  0  1 

“0  •  •  •  0" 

^>Vi 

and  B  = 

o'. 

A  ^  A 

0  .  .  .  K 

Now  we  can  derive  equations  for  x^  in  terms  of  x^_£.  From  (4.1)  we  get 

■1  =  A^ 

(A10) 


Xm  ,•  =  Ax..  .  ,  +  BA,.  . 


(A12) 


=  Ax„ 


-i-i  ♦  »<  i0(AT)1'Jfi-j  *  (aT)Ww 


‘"V,,  +  Bj0(#T)i'JcV1(zN-3  -  CW  *  ^^Vl 


For  convenience  define 
(A13) 


v,  -  [B(ATP]m  =  KCAT]; 


-  A7  - 


and 

(A14)  Wi  =  ViCTR_1  , 


where  the  subscript  m  indicates  the  block  of  the  last  m  rows. 

Note  that  is  a  matrix  of  order  mxn  which  can  be  computed  recursively. 
Taking  advantage  of  the  structure  of  A  and  B  no  more  than  n*m  operations  are 
necessary  for  each  V... 

Thus  using  (A13),  { AT 4 )  and  (4.6),  equation  (A12)  gives  the  relation  for 
the  last  m  components  x^_ .  of  x^_^  as 


<»'5>  XN-1  ■ 


4XN-j-y+4  +  Vi+1^N+1  • 


Before  collecting  terms,  let  k  =  -j-y,  then 


xN-i  =  J^Vn-1-I-y+I  +  j^0Wi-jZN-j  ' 


(A16) 


1  f-, 

<=-l-Y  4=1  ' 


+k+yC4XN+k+4  +  Vi+l^i+l  * 


Now  apply  the  index  transformation  (AT )  with  =  1 ,  j'2  =  y,  =  -i-y, 
i2  =  -Y-  This  yields  an  expression  of  the  form  (k  =  j-4) 


-Y  i+j+y  -i  j+i+y  0  Y 

{A1?)  k  4  j=l -i-y  4=1  aj"£’Sl  +  j  =  5-y  4=^“***  +  j  =  -Hl 


Also,  since  A  =  -lm  we  have 


(A18) 


-x 


N-i 


4=1 


A£,XN-i  -1  -y+4 


Y+l_ 

=  l  A„x 

4=1 


4  N-i-l-y+4 


If  we  let  j  =  4-i-l  this  becomes 


(A19) 


Y+l  Y-i 

J^4XN-i-Y+4-l  =  jJ_iRj+1+lXN-Y+j  ‘ 


Now  we  can  use  (A17)  and  (A19)  to  rewrite  (A16): 


Y-i 


°  jJ_i  AJ+i+'>XN-Y+j  +  Wi-jZN-j 


b-f-1  il  'WtCtWj 


.  + 


(A20) 


+  %  Jjj  Wi’Wj  + 


+  jf  j  U  cx 

j=y-Zi+1  i*j  ‘+j“^^^XN"Y+j 


+  Vi+1^N+1 


This  holds  for  0  <  i  <  Y-l.  After  collecting  terms  involving  for 


j  >  1  on  the  left  hand  side,  we  obtain  a  system  of  equations  of  the  form 
(A21 )  tV  =  T1^  +  zl 


where 


= 


Vy+1 


»  X 


N-2y+l 


'N-y 


i —  ct  —i 


,  zl- 


and 


(A22 ) 


Zi+1  "  Wi-jZN-j  +  Vi+l^l+l  »  1  0  •••  Y-l 


and  the  matrix  T®  is  given  by 


(A23) 


T°  .  = 
i.J 


j+i-1 


"Ai+j  +  Jj  Wi-l+j-R.v'K. 


l  4=j 


U 


1  +  j  -  JlC  ft.  ’ 


i_ec#  .  If  1  «  j  S  Y-i+1 

if  y-i+2  S  j  s  Y 


This  can  be  defined  recursively,  too.  Let 


-  A9  - 


f  i+j-1 


(A24 ) 


FU 


1  Wi+j-£-l  CJl  ’  1  “  J  -  Y" 


*=j 


i+1 


j.  wi-i+j-A  ’  Y'i+2 


Then  for  1  <  j  <  y- i+1 


Fi+1  J-1  =  Wi  +J-JL-1  C£  =  Wi C  j  - 1  +  Fi,j 


Similarly  for  j  >  y-i+1 : 


Fi+1 


.J'1  #JL  Wi- 


£=j_l  "i-l+j-A  =  WiCj-l  +  Fi,j  ’ 


Consequently  F  satisfies  the  recursion 


(A25) 


Fi+l,j-l  Fi , j  +  WiCj-l 


in  both  cases.  Thus,  after  F.  .  and  F.  have  been  computed,  the  recursion  is 

l  i >Y 

initialized,  and 

(A26) 


T0  f  -Ai+j  +  Fi,j  if  1  *  J  *  ^i+1 
1  ,J  F.  .  otherwise 

'  9  J 


rl 


To  find  T  change  the  index  on  the  right  hand  side  from  j  to  k-y.  Then  we 
get 

J,  ^i+1  -y+kXN-2y+k  '  k=Ji+y  Wi+k-y- AXN-2y+k  ’  0  -  1  ‘  Y_1 


From  this  we  see  that  for  1  i  i  s  y 


(  0 
A 


(A27)  Ti  ,|, 


if  k  <  y-i+1 
if  k  =  y-i+1 


k-y+i-1 


A1-y*k  -  T,  “i-Uk-y-A  if  k  =•  Y-(*1 


-  A10  - 


For  k  >  Y_i+1  define 


_  k-Y+i-1  _ 

Ti.k  =  ^  Wi-l+k-Y-«.C«- 


Then  clearly 
(A28) 


Ti,k  “  Ti-1  ,k+l  ’ 


and  we  only  need  to  compute  T.  for  i  £  2,  i.e., 

*  >Y 


i-1 


(A29) 


i.Y  =  J, 


Z=1 


Z  ’ 


APPENDIX  B 


INVERSION  OF  BLOCK  CIRCULANT  MATRICES 
A  matrix  of  the  form 


(Bl)  Hc  - 

where  itself  is  a  matrix  of  dimension  dxd,  is  called 
matrix  with  dimension  N  and  base  dimension  d. 

Consider  the  first  column  of  blocks  HQ,  Hn_-j  »  ...  H 
of  length  N  by  taking  the  elements  in  the  position  (i ,j) 
Nq  (i»j)»  H  -j  ( i » J  ) »  •••  H-j  (i  ,j ) ). 

Then  the  inverse  of  Hc  can  be  computed  efficiently 
Step  1 :  Take  the  DFT  of  the  sequences  Hg(i,j),  Hn_-j  ( i ,  j 
(i,j)  and  denote  those  transforms  by  Hq ( i » j ) ,  HN_-|(i»j), 
Step  2:  Compute  the  inverses 

.  .  .-1 

Bk  =  Hk  ,  for  k  =  0,  ...  N-l 

Step  3;  Take  the  inverse  DFT  of  the  sequences 

Bg(  i  i  j  )  i  *'•  » J  ) 

and  denote  the  result  by 


a  block  circulant 

,  and  form  sequences 
of  each  matrix,  i .e. , 

in  the  following  manner: 
...  H-|  ( i , j )  for  each 
...  H^(i ,j). 


B2  - 


which  gives  the  first  column  of  blocks  of  H~^. 

c 

Step  4:  Circulate  the  blocks  Bk  to  obtain  the  matrix  B  =  H”1  given  by 


Proof:  Since  H  is  a  NxN  block  circulant  matrix  with  base  d,  H  can  be 
^  c 

expressed  as  a  sum  of  Kronecker  product  of  matrices  as 


(B3) 


(B4)  CQ  = 


N-l 

Hc  »  l  ck©Hk 

c  k=0  K  K 


»  C-,  = 


'0.1 


>  »  •  •  ,  C 


N-l 


0  1 


1  0 


where  each  Ck  is  NxN  and  each  is  dxd.  Define 


Hc  =  F  Hc  F* 


(B5) 


N-l 


‘  (F  ©  Id)l  l  Ck  ©  Hk)(F*  ©  Id) 
k-0 


where  F  is  the  NxN  unitary  discrete  Fourier  transform  matrix.  Then  by  properties 
of  the  cross  project  of  matrices,  we  get 


(B6) 

(87) 

(B8) 


N-l 

Jo 


N-l 

Jo 


N-l 

l 

k=0 


(F  ©  Id)(Ck  ©  Hk)(F*  ©  Id) 

(F  Ck  ©  Hk)(F*  ©  Id) 

(F  Ck  F*)  ©  Hk 


-  63  - 


Since  Ck  is  a  circulant  matrix  with  l's  and  zeros,  the  product  F  C1  F* 
is  simply  a  diagonal  matrix,  Dk>  i.e.. 


(B9) 


Dk  =  F  Ck 


dk(0)) 


dk(D 


0 


0 


dk(N-l) 


-i  —  ke 

where  dk(£)  =  e  J  N  for  k,£  =  0,1,  ...  N-l 

(BIO) 

N-l 

Hence  H  =  £  D.  <x)  H,  (Bll) 

c  k=0  K  K 

A  A 

Hc  is  a  block  diagonal  matrix  with  the  diagonal  blocks  as  H^,  i.e., 


A 


N-l 


<V*>Hk  • 


Substituting  (BIO)  results  in 

N-l  ,  ?!  k 

H  (i,j)  =  l  e  3  N  K\(i,j) 

*  k=0  K 

A  A 

Since  Hc  is  a  block  diagonal  matrix,  the  inverse  of  Hc  is  simply  the  inverse 
of  each  block  matrix  on  the  main  diagonal. 

In  order  to  get  B,  the  inverse  of  Hc,  the  inverse  DFT  is  required. 

Note  that 

(He)"1  =  (F  Hc  F*)"1 
"c  =  (F*)~1HC_1 (F)_1 

H;1  =  F*H~]  F  =  F"1H"1(F"1)* 


^  -  •  - 


The  inverse  of  a  block  circulant  matrix  is  also  a  block  circulant  matrix. 


This  can  be  shown  as  follows: 


8  can  be  expressed  as  §  =  £  E.  ©  B. 

k=0  K  K 


where 

"1 

0 

0 

ro  o 

0  1 

0 

11 

O 

UJ 

0 

— 1 

•  E1 = 

:  o 
% 

4  •••’  EN-1 

°  •* 

—1  o 

1 _ 

1  N-l 

Then  B  =  H'1  =  (F*  ©  Id)(  l  Ek  ©  Bk)(F  ©  Id) 

k-0 


N-l  _  - 

=  l  (F*  E.  F)  ©  B 
k=0  K  K 

It  is  found  that  F*  F  is  a  circulant  matrix.  Therefore  (F*  E^F)  @  B 
is  a  block  circulant  matrix.  Because  summation  of  block  circulant  matrices 
gives  a  block  circulant;  B  is  a  block  circulant  matrix. 
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